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Abstract 

By subjecting a dynamical system to a series of short pulses and varying several time delays we 
can obtain multidimensional characteristic measures of the system. Multidimensional Kullback- 
Leibler response function (KLRF), which are based on the Kullback-Leibler distance between the 
initial and final states, are defined. We compare the KLRF, which are nonlinear in the prob- 
ability density, with ordinary response functions (ORF) obtained from the expectation value of 
a dynamical variable, which are linear. We show that the KLRF encode different level of infor- 
mation regarding the system's dynamics. For overdamped stochastic dynamics two dimensional 
KLRF shows a qualitatively different variation with the time delays between pulses, depending on 
whether the system is initially in a steady state, or in thermal equilibrium. 

PACS numbers: 05.40.-a, 89.70.Cf, 02.50.Ey 
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I. INTRODUCTION 



One of the most common ways to investigate the properties of a dynamical system is to 
study how it responds to controlled external perturbations. The response of a system to a 
weak perturbing field is related to its equilibrium fluctuations by the celebrated fluctuation 
dissipation relation jl, 2|. The response provides a direct measure of system dynamics and 
fluctuations. 

In a time-domain response measurement one uses a series of impulsive perturbations 
(Fig. [1]) and records some property of the system as a function of their time-delays. Impulsive 
perturbations make it possible to study the free dynamical evolution of the system during 
the time delays unmasked by the time proflle of the perturbing field. Furthermore, the joint 
dependence on several time delays can be used to separate the contributions of different 
dynamical pathways. Due to its dependence on multiple time delays this method is termed 
mult idimensional . 

The response of a system is typically measured by the expectation value of some op- 
erator l|-l3|. This is a linear functional of the system probability density (or the Density 
matrix for quantum systems). Multidimensional response have had considerable success in 
nonlinear spectroscopy, due to the ability to control and shape optical fields. Applications 
range from spin dynamics in NMR |4j, vibrational dynamics of proteins in infrared systems 
and electronic energy transfer in photosynthetic complexes as probed by visible pulses 5|, |6|. 
These span a broad range of timescales from milliseconds to femtoseconds. 

Interestingly, there exist non-linear functionals of probability densities which have in- 
teresting physical interpretations. One such quantity is the Von-Newman entropy S{p) = 
— Tr p In p. A related quantum nonlinear measure called the concurrence serves as a measure 

nn 

of quantum entanglement [7|, |8|. The Kullback-Leibler distance (KLD) or relative entropy, 
Tr po In Po/P) which compares one probability distribution to another, is a nonlinear measure 
that had been found useful in many applications. This paper aims at developing multidi- 
mensional measures based on the KLD. 

Numerous applications of the KLD differ in the probability distributions involved. The 
ratio of the probability of a stochastic path and its reverse at a steady state has been 
connected to a change of entropy j9-12|. For an externally driven system a similar quantity 
was found to be related to the work done on the system [131]. As a result, the KLD which 
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compares the path distribution to a distribution of reversed paths is a measure of the lack 
of reversibihty of a thermodynamical process. 

For distributions in phase-space (as opposed to path-space) j^_the KLD between the density 



of a driven system and the density of a reversed process 
the driven density and the corresponding equihbrium density 



, or the distance between 
15| (for the same value of 



parameters) were shown to be bounded by dissipated work in the process. The transfer of 
information through a stochastic resonance is quantified by the KLD between the probability 
distributions with and without the external input The ability of neuronal networks 



to retain information about past events was characterized by the Fisher- information |18| . 
which is closely related to the KLD between a distribution and the one obtained from it by 
a small perturbation. 

We shall examine the response of a system to impulsive perturbations which drive it 
out of a stationary (steady state or equilibrium) state. The KLD between the distribution 
before and after the perturbation does not correspond to an entropy, or work. However, 
since it compares the perturbed and unperturbed densities, it characterize how "easy" it 
is to drive the system away from its initial state. In ordinary response theory, one com- 
pared the expectation values of some operator taken over the perturbed and unperturbed 
probability density. This depends on the specific properties of the observed operator. The 
KLD is a more robust measure for the effect of the impulsive perturbation on the probabil- 
ity density. By expanding the KLD in the perturbation strength we obtain a hierarchy of 
Kullback-Leibler response functions (KLRF). These differ qualitatively from the hierarchy 
of ordinary response functions (ORF), since they are nonlinear in the probability density. 
The KLRF serve as a new type of measures characterizing the dynamics and encoding dif- 
ferent information than the ORF. For example, the second order KLRF, which we connect 
to the Fisher-information, is found to exhibit qualitatively different dependence on the time 
delays, depending on whether the system is perturbed out of a steady state or out of ther- 
mal equilibrium. This is in contrast to the corresponding ORF. The fluctuation dissipation 
relation jj, |2| , which is linear in the density matrix, can also distinguish between systems 
driven out-of-equilibrium and out of a steady state. The KLD offer a different window into 
this aspect. 

The structure of the paper is as follows. In Sec. [TTl we describe the multidimensional 
measures and present the two heirarchies of ORF and KLRF response functions. These are 
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then calculated using a formal perturbation theory in the coupling strength to the external 
perturbation. In Sec. IIIII we show that for systems undergoing overdamped stochastic dy- 
namics the non-linear KLRF are naturally described using a combination of the stochastic 
dynamics and its dual dynamics. In Sec. HV] we extend the results of Sec. IIIII to discrete 
Markovian systems with a finite number of states. Our results are discussed in Sec. El 

II. MULTIDIMENSIONAL MEASURES FOR NONLINEAR RESPONSE BASED 
ON THE KULLBACK-LEIBLER DISTANCE 

We consider a system initially at a stationary state (either equilibrium or a steady state), 
which is perturbed by a series of short pulses, as depicted in Fig. [1] The probability dis- 

Si S2 




FIG. 1: Heuristic schematic of the process studied. A system is prepared at an initial probability 
distribution po and subjected to impulsive perturbations. The kth pulse is centered at and its 
strength is denoted by s^- tj = Tj+i — tj are the time intervals between successive pulses. 

tribution describing the driven system at time r, p(t), depends parametrically on Si, the 
strength of the ith pulse, as well as on the time differences between the pulses ti = Xj+i — Tj 

with Tn+l = T. 

Ordinary response theory focuses on the expectation value of some observable 

{A) (t) = Tr [p{t)A] , (1) 
and its dependence on the parameters Sj, r^. The lowest ordinary response functions (ORF) 
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and so forth. The time differences in Eqs. (E])-® can be expanded in terms of the time 
delays between the pulses, r — = J2k=i^k- T^^''^ are used to investigate various properties 
of the unperturbed dynamics, such as the existence of excited modes, and the relaxation 
back to a steady state. 

Here, we focus on different, but closely related quantity. Instead of studying an expec- 
tation value of an observable, we focus on a quantity that compares the perturbed and 
unperturbed probability distributions. The KLD, also known as the relative entropy, is 
defined as 



^^(PoIIp) = Trpoln 



Po 
P 



ln^\ 
P /o 



(5) 



The KLD vanishes when the two distributions are equal V {p\\p) = 0, and is positive oth- 
erwise, V (p'Wp) > for p'{x) 7^ p{x) fl9|. Note that the KLD is not a true distance since 
^ (p'IIp) 7^ ^ (pIIpO it does not satisfy the triangle inequality. 

The KLD measures the dissimilarity between two distributions. It had found many appli- 
cations in the field of information theory 19|]. For instance, the mutual information between 
two random variables x,y is V {P{x,y)\\P{x)P{y)), where P{x,y) is the joint distribution 
while P{x) {P{y)) is the marginal distribution of x (y). In the present application the 
Kullback-Leibler distance is a measure for the deviation of the system from its initial state. 

In a manner similar to the definition of the ORF, we define a KLRF hierarchy by taking 
derivatives of the KLD with respect to pulse strengths, and displaying them with respect to 
the time delays 
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All the derivatives are calculated at s = 0, and we have used the relation lims-i>o p{t) = Po- 
This is also true for all other s derivatives in the following. To keep the notation simple we 
will not state this explicitly. Higher order KLRF are defined similarly. It is important to 
note that the ORF are linear in p whereas the KLD are nonlinear. We thus expect the KLD 
to carry qualitatively different information about the dynamics. 

The second derivative ([7]), known as the Fisher Memory (or information) matrix, plays 
an important role in information theory, since the Cramer-'Rao inequality means that it is a 



measure of the minimum error in estimating the value of a parameter of a distribution |19 |. 
The Fisher information has been used recently to analyze the survival of information in 
stochastic networks fl8|. 

Conservation of probability implies that the first KLRF vanishes 

Q^\r -Ti) = - J dxpo-^ lnp = - J dx-^p = 0. (9) 

The second derivative, the Fisher memory matrix, is given by 

^(2), s /", _idp dp / d\npd\np \ 

A straightforward calculation allows to recast the third order derivative of V in terms of 
products of lower order derivatives 

-ThTi- Tj, Tj - Tfc) = 

dlnp d'^\np\ ^ / din p d'^ln p\ ^ / din p d'^\io.p\ ^ / din p din p din p\ ^^^^ 
dsi dsjdsk/Q \ dsj dsidsk/^ \ dsk dsidsj / ^ \ ds. dsj dsk / q' 

In what follows the derivatives will be calculated perturbatively in sj. It is important to 
note that the A^'th derivative of In p has contributions from interaction with at most 
pulses. The contribution from the linear component, which interacts with A^ pulses, has the 
same structure of the perturbation theory for observables (which is also linear). However, 
since In p is a non linear function of p the A^'th derivative contains a non-linear contribution 
which is a product of lower order contributions for p. The KLRF encode qualitatively new 
information about the system dynamics in comparison to the ORF. 
The time evolution of the probability distribution is given by 

^ = -^(-)P- (12) 
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This formal equation is quite general, and can describe either Hamiltonian (Unitary) or 
stochastic dynamics where the operator £ will accordingly be the Liouville, or the Fokker- 
Planck operator. 

For a system subjected to a time dependent weak perturbation we can write 

t{T)=t, + t'{T), (13) 

where we assume that the unperturbed system is time independent and is intially in a steady 
state, so that JCqPq = 0. We consider an impulsive perturbation of the form 

n 

1=1 

where Ca describes the action of a pulse on the probability distribution, and Sj is the overall 
strength of the ith pulse. 

Using these definitions, the state of the system at time r can be expanded as a power 
series in the number of interactions with the pulses 

p(t) = po + 5pi + Sp2^ (15) 

i=l i=l j=l ^ ^ 

The partial corrections for the density, S^^\ appearing in Eq. fllSp . contain all the informa- 
tion necessary for computing both the KLRF and ORF. They are given by 

S^-^\x;t) = J dxidx2l(o{x,X2;t)CA{x2,Xi)po{xi), (16) 
S^'^\x;t',t) = j dxi - ■ ■dx4Jo{x,Xi]t')CA{xi,x-i)Uo{x-i,X2\t)CA{x2,xi)po{xi). (17) 

Here Uoij — r') = exp — (r — t')Cq is the free propagator of the unperturbed system. 
Conservation of probabihty requires that / dxSpi{x), which, in turn, means that J dxS^^^ = 
0. 

Eqs. ffT^ - ffT7|) can be used to calculate the logarithmic derivatives, which then determine 
the KLRFs. We will only need the first two logarithmic derivatives, which are given by 



d\np 1 d5pi 
dsi po dsi po{x) 



^5«(x;r-7:,), (18) 
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and 

dlnp 1 d'^6p2 1 d6pid6pi 



dsidsj po dsidsj p^ dsi ds 



^-5(2) (x; r - n, n - Tj) - -^S^^\x- r - n)S^^\x- r - r,). (19) 
Pq{x) pg(x) 



To calculate the ORF, we substitute Eq. ( 1151) in Eq. ([T]), resulting in 

{A}^ = (A)o + ^ / dTn / dTn-1 ■ ■ dTiS{Tn) ' " " s{Ti)71^'^\t -Tn, Tn-Tn-1, " " " , T2-n). 

n=l "^"^0 J TO J TO 

(20) 

It is interesting to compare the expressions of the KLRF the ORF. We Calculate Q^^^ 
and 71^^^ for j perturbing pulses. To leading order, we find 

Q?ih) = 0, 

7^(^)(tl) = j dxA{x)S^^\x;ti). (21) 
At the next order, we compare the Fisher information to the second order ORF, 

Quit2,ti) = j dxpo\x)S^'\x;h+t2)S^'\x;t2), (22) 

and 

7^(2)(t2,^l) = j dxA{x)S^^\x;t2,ti). (23) 

(2) 

The non-diagonal elements of Q]j depend on the two delay times. Expressions for the third 
order response functions are given in App. IB 1[ 

Both TZ^^"^ and Q*--'-' depend on the same set of j time intervals with some important 
differences. Q^^'^ vanishes, while the linear response T^*^^-* does not. Q^'^^ and TZ^'^^ have a 
different structure: TZ^'^^ can be calculated from the second order correction to the density (or 
S^"^^) while Q'-^^ is determined from a product of iS'^^^'s describing the first order interaction 
with different pulses. This difference reflects the non-linear dependence of the KLRF on p, 
and also applies to higher orders. 

A comment is now in order regarding our choice of the KLD (JSj). We have chosen to 
use V (poIIp) as the measure for the effect of the perturbations. V (p||po) would have been 
equally suitable. However, as discussed in App. El the leading order of both KLDs in the 
strength of the perturbation, i.e. their Fisher informations, coincide. Therefore all the 
following results pertaining to the Fisher information would hold for either choice. 
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III. APPLICATION TO OVERDAMPED STOCHASTIC DYNAMICS 



In the following we use the formal results of Sec. [IT] to calculate the leading order ORF 
and KLRF for a system undergoing overdamped stochastic dynamics. We show that the 
Fisher information is related to a forward-backward stochastic process. The backward part 
is driven by the po-dual process, which will be simply referred to as the dual in what follows. 
The Fisher information is found to exhibit qualitatively different properties for systems 
perturbed from equilibrium, or from a steady state. We also use the eigenfunctions and 
eigenvalues of the dynamics to derive explicit expressions for several low order ORF and 
KLRF. 

In stochastic dynamics the probability density plays the role of a reduced density matrix, 
which depends on a few collective coordinates. In this reduced description the entropy 
— Tr plnp typically increases with time. This should be contrasted with a description which 
includes all the degrees of freedom, where the dynamics is unitary and the entropy does not 
change in time. For completeness unitary dynamics is discussed in App. |X1 



A. Fisher information for systems perturbed out of equilibrium vs steady-state 

The Fisher information can be represented in terms of the dual stochastic dynamics. 
This interesting property reflects its non-linear dependence on p. We examine a stochastic 
dynamics of several variables x^, given by 

| = i^.(x) + e.(t;x). (24) 

Here we use the Ito stochastic calculus. ^ The noise terms are assumed to be Gaussian with 

Ut)Ut')) = G^At-t')- (25) 

with G'(x) a symmetric positive definite matrix. While for many systems this matrix does 
not depend on the coordinate, Gij = ^ksTdij, this assumption will not be used in what 
follows. 



^ The Ito and Stratonovich calculus offer two different recipes of interpreting Eq. (f24| . Both methods are 



equally viable as long as they are used in a consistent manner. Details can be found is Ref. 22 1 
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Equation is equivalent to the Fokker-Planck equation 

|--E7^^<P + Ei7dLG,p^-£„p. (26) 



dt ^ dxi ^ 2 dxidx. 
t ij 



In what follows we present the dual dynamics, which can be loosely thought as the time 
reversed dynamics: it have the same steady state, but with reversed steady state current. 
We consider the current density 

J^ = F^P-Y.\-i^.G^^P- (27) 



The Fokker-Planck equation can be written in terms of the current 

dp _ ^ dJj 



dp dJi 

Y.7dr.- (28) 



The steady-state is the solution of 

t^Po = 0. (29) 

We write 

Po = e-^°, (30) 

which defines 0o- For systems at equilibrium 0o is simply the potential. However, this is 
not the case for general steady states. The steady state current can be written as 

4"' - ^'P^ - E 4 E G« (E ^OJ^^^ - E Oji^ .0. (31) 

j ■' j \ k kl / 

After some algebra, the generator of the stochastic dynamics can be written in terms of the 
steady state density and currents {20! 

^o^-E|-^n'°V-iE|-««-*^^-p- (32) 

The dual dynamics is given by 
with 

t'' = potlp,\ (34) 

A straightforward calculation gives 

^'p-El-^V-EiaSF'^.p (35) 

10 



with 

Ft = -F. + Y.et>^^G,,e-^^. (36) 

It is a simple matter to verify that the dual dynamics has the same steady state as the 
original one, but the steady state currents have opposite signs. It can be simulated by 
integrating the Ito stochastic equation 

^ = Ft{^) + Ut;^). (37) 

The dual dynamics reverses the non conservative forces in the system. This relates the joint 
probability to go from one place to another in the original dynamics to the joint probability 
of the reversed sequence of events in the dual dynamics 2^ 

Po(x',ti|x)po(x) = P^(x,ti|x>o(x'). (38) 

The left hand side of Eq. ( l38i) is the joint steady state probability to first the system at 
X, and at x' after a time ti. The right hand side is the joint probability of the reversed 
sequence of events, but for a modified dynamics. When this modified dynamics is the dual 
these joint probabilities become equal. 

We next turn to discuss the system's response to a series of impulsive perturbation. We 
assume that the perturbation is of the form 

£^p=lyAMp (39) 

7 ^ dxi dxi 

with A{x) as a potential field perturbing the system. Using Eq. f ll6p we have 

5'i(xi;ti)= /"cixoPo(xi,ti|xo)-5(xo)po(xo), (40) 
J 7 

where 

^W = i:^-^^- (41) 
With the help of equation ( 138|) . we obtain 

5i(xi,ti)= / dxo-P(xo)P'^(xo,ti|xi)po(xi). (42) 
J 7 

We now have all the tools needed to compare the ORE and KLRE for overdamped 
stochastic dynamics. The leading order response function is given by 

^i(^i) = - / (ixio?xoyl(xi)Po(xi,ti|xo)P(xo)po(xo), (43) 
7 J 
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while Q^^^ = 0. At the next order, we have 

Qi2ih,ti) = ^ J c?X2rfxirfxo5(x2)P''(x2,t2|xi)Po(xi,ti + ^2 |xo)5(xo)po (xq) , (44) 



and 

1 f d { OA 

Tl2{t2,tl) = J '^X2(iXi(iXoA(x2)Po(x2,t2|xi)^^ ^ 



Po(xi,ti|xo) 



P(xo)po(xo). 

(45) 

Some insight into the structure of different response functions can be gained by represent- 
ing them as ensemble averages over stochastic trajectories. The first order response function 
can be simulated directly using stochastic trajectories of the original dynamics. The ap- 
pearance of a derivative of the conditional probability complicates the direct simulation of 
7^2- It may be possible to circumvent this difficulty using the finite field method, where one 
combined simulations with and without a finite, but small perturbation [21j. can be 
simulated with trajectories which follow the original dynamics for time ti + ^2 and then the 
dual dynamics for time ^2- 

Systems at equilibrium are self dual, allowing to substitute / (ixiP'^(x2, t2|xi)Po(xi, ti + 
^2|xo)) = Po(x2,ti + 2t2|xo)in Eq. (I44l) . As a result the Fisher information only depends 
on a single time variable ti + 2t2- This is in contrast to systems which are perturbed from 
a nonequilibrium steady state, whose Fisher information is a two dimensional function of 
ti and ^2- The Fisher information is therefore qualitatively different for systems which are 
perturbed out of a steady state, or out of an equilibrium state. 

For the self-dual case Q^^2 ^^.s the same structure as TZ''^\ up to a replacement of A{x) 

(2) 

with B{x). However, in the general, non self-dual case, the structure of is manifestly 
different than that of 7^(^) and TZ^'^\ 



B. Eigenfunction expansion of the Fisher information 

An alternative approach for the calculation of the Fisher information, as well as higher 
order KLRF, uses eigenfunction expansions for the density p. It will be sufficient to examine 
a simple one dimensional model where Lq is a Fokker-Planck operator and the perturbation 
is given by Eq. fl39l) . 

The propagation of the unperturbed system can be described in terms of the eigenvalues 
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and eigenf unctions of Cq. The right eigenfunctions satisfy 

CqPu = J^{n)pn- (46) 

Similarly, the left eigenfunctions, g„, satisfy g„£o = ^{n)qn. It is assumed that the right 
and left eigenfunctions constitute a byorthogonal system, that is 



dxqn{x)pm{x) = 6nm- (47) 

Any probability density can be expanded in terms of right eigenfunctions, 

p{x) = J2CnPn{x), (48) 

n 

with 

Cn = y dxqn{x)p{x). (49) 

We consider systems perturbed out of thermal equilibrium, with a probability density 
Pq{x). In this case the left eigenvalues are simply related to the right eigenvalues. For 



variable x which is even with respect to time reversal, this relation takes the form 



22 



Pn{x) = qn{x)po{x). (50) 

As a result, one can use only the left eigenfunctions, which in this case satisfy the orthogo- 
nality condition 

j dxqn{x)qm{x)po{x) = 5nm- (51) 

We note that qo{x) = 1 and A(0) = correspond to the equilibrium distribution. 

Our goal is to calculate the response of the system to a perturbation with a coordinate 
dependent operator A{x). Several types of integrals appear repeatedly in the calculation, 
and it will be convenient to introduce an appropriate notation. One comes from the need to 
decompose the probability distribution into eigenstates after each interaction with a pulse 

/d IdA 1 f dqm dA 

dxqm{x)— -^qn{x)po{x) = - dx-^—qn{x)po{x), (52) 

where integration by parts was used in the second equality. (We assume that po{x) falls of 
fast enough to eliminate boundary terms.) The calculation of the response also involves an 
evaluation of the average of an observable, which, in the current setting, leads to integrals 
of the form 

Cau = J dxA{x)qn{x)po{x). (53) 
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It is straightforward to calculate S^^^ and S^'^^ with the help of Eq. (I52p . We find 

5i(x;ti) = ^S„oe-^(")*^g„(x)po(a;), (54) 

n 

and 

5(2)(x;t2,ti) = ^S^„S„oe-^(")*^-^("^)*^g^(x)po(x). (55) 

mn 

We now calculate the first few ORF and KLRF. The first order response functions is 

7^l(^l) = 5^CA„i3„oe-^(")*^ (56) 

n 

Similarly, the second order response function is given by 

7^2(^2,^l) = J^c^^s^.is^oe^^^")*^-^^'")*^ (57) 

mn 

TiS^'' should be compared to the KLRF of the same order, namely the Fisher information 

Qg(t2,ti) = $^S^oe-^(")[*^+^*^l, (58) 

n 

which is calculated with the help of Eqs. fl22|) and fl5T]) . Again, the Fisher information of 
systems perturbed out of equilibrium depends only on ti + 2t2. 
For a system initially in a steady state the Fisher information is 

nm 

However, since the relation Eq. (!50|) does not hold in this case, the integral in Eq. (!59|) does 
not vanish for n m, and the Fisher information becomes a two dimensional function of ti 
and t2. 

are calculated in App. lB 2l where it is shown that Q^i23 is a three dimensional 
function of all its time variables. This qualitatively different signature of systems initially 
at steady state vs equilibrium is unique to the Fisher information, due to its quadratic 
dependence on 6 p. It does not apply to higher order quantities, such as Q^^\ 

The eigenf unctions for a simple example, of an harmonic oscillator with an exponential 
perturbation, are presented in App. O 
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IV. MASTER EQUATIONS, DUAL DYNAMICS, AND THE FISHER INFORMA- 
TION FOR DISCRETE SYSTEMS 



The description of the KLRF in terms of a combination of the regular stochastic dynamics 
and its dual holds also for Markovian systems with a finite number of states. Below we derive 
a simple expression for the Fisher information of a stochastic jump process in terms of the 
dual dynamics of the original process. 

A. Dual dynamics of discrete Markovian systems 

Consider a system with a finite number of states, undergoing a Markovian stochastic 
jump process described by the master equation 



P is a vector of probabilities to find the system in its states and TZ is the transition rate 
matrix Its off diagonal elements are positive, TZij > for z 7^ j, and express the rate of 
transitions from state j to i, given that the system is at j. The diagonal elements satisfy 
J2i'^n — 0- further assume that there exist a unique steady state, P^'*^ satisfying 



and that at this steady state there is a non vanishing probability to find the system at each 
of the states i, P/*'' 7^ 0. The master equation is one of the simplest models for irreversible 
stochastic dynamics. Below we briefly describe some of its relevant properties, such as the 
backward equation, and its dual dynamics. 
One can define an evolution operator 



p = 7^p. 



(60) 



7^p(") = 



(61) 



(62) 



which satisfy the equation of motion 




(63) 



with the initial condition Ui 



■0 = 



1. 
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For this model a dynamical variable A is a vector with a value corresponding to each 
state of the system. Its expectation value is given by 

(A),^A-P(t) = 5;]APKi). (64) 

i 

Instead of calculating this average by propagating P, one can define a time dependent 
dynamical variables, A(t), such that (A)^ = A{t) ■ P(0). This is the analogue of the 
Scrodinger and Heisenberg pictures in Quantum Mechanics. The equation of motion for 
A(i) is easily shown to be 

|A(t) = 7^tA(^). (65) 

This equation is known as the backward equation, which turns to be related to the dual 
dynamics. (The backward equation is also written as a function of an initial rather than a 
final time.) TZ and V) have the same eigenvalues. However, the roles of the right and left 

eigenvectors are interchanged and A(t) decays to a uniform vector with at long times. 

is) 

Let us now define a diagonal matrix 11, using liu = . The dual evolution is then 
defined as 

= ^7^^^-^ (66) 

The fact that IT is built from the steady state of 7^, P'^, guaranties that TZ*^ is a physically 
reasonable rate matrix, that is, that is satisfies Yli'^ij — 0- The dual dynamics describes 
a physically allowed process which has the same steady state as the original process it was 
derived from. However, at this steady state the dual currents have opposite signs compared 
to the steady state currents of the original dynamics. A process is self-dual, that is, Tif" — TZ, 
if and only if it satisfies detailed balance. Self duahty is therefore related to being in thermal 
equilibrium. 

B. Perturbation theory and the Fisher information 

Consider a system subjected to several impulsive perturbations 

n{T) ='R^ + Y^ sJIa5 (r - Ta) . (67) 

a 

Where TZa corresponds to some physical perturbation A. Here the free evolution is given by 

= e^°*. (68) 
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Similarly the evolution during an impulsive perturbation can be described by 

P(r, + e) = e^"^^P(r„-e), (69) 

where e is arbitrarily small. 

We expand the exponent in Eq. ( !69l) . and collect all terms of the same order in Sa- The 
expression for the Fisher information, following Eq. ( 122|) is 

Q?^{t2,h) = Y.^K'^aP')^ {ul^,^nAP')^. (70) 

i i 

This expression can be simplified by writing one of the propagators in terms of the dual 
process, using the relation 

g7^^^^-l = n-^e^'*. (71) 
After some algebraic manipulations, we find 

Qf^it^M) = (7^1e^'*^e^«(*^+*^)7^^P^)^. (72) 

i 

In Eq. (ff2l) we have defined TZj^ = IIT^^^n^^ TZj^ is not the dual of TZa since 11 is not 
composed of the eigenvector of TZa which corresponds to a vanishing eigenvalues. (11 is 
composed of the eigenvalue of T^o-) 

Equation (172|) is the analogue of Eq. (l44l) for discrete systems. It demonstrates that 
for this model the fisher information can be simulated using a combination of the ordinary 
process and its dual. However one must also include a (possibly artificial) dual perturbation, 
TZj^, which can nevertheless be computed from the known physical one. As before. Equation 
(1721) shows that the Fisher information, for self-dual systems, is a function of ti + 2^2- 

As a side remark, the model considered here has only even 'degrees of freedom' under time 
reversal. More general models may also include odd degrees of freedom, such as momenta. 
In that case one may speculate that there would be anti-self-dual systems whose dual turns 
out to be the time reversed dynamics. This would lead to a Fisher information depending 
on ti alone, as is the case for Unitary dynamics. 

V. DISCUSSION 

In this work we have studied a system driven from its steady state by a sequence of 
impulsive perturbations. We have defined a new set of measures for the response of the 
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system to the perturbation, the KLRF, which are given by the series expansion of the 
Kullback-Leibler distance between the perturbed and unperturbed probabihty distributions. 

At each order the KLRF and ORE depend on the same time differences between the 
pulses. However, there are important differences stemming from the nonhnear dependence 
of the KLRF on p. The expression for the KLRF, for instance Eqs. (12^ and (]B2p reveal 
quantities which can be simulated using several trajectories which end at the same point. We 
have shown that a simpler, but equivalent description exists. It uses the dual dynamics which 
allows to "run some of the trajectories backward". This description is especially appealing 
for the Fisher information. Instead of viewing the Fisher information as composed over sum 
of pair of trajectories joined at their end point one can view it as an average of contributions 
of a single forward-backward trajectory. 

Another difference between the KLRF and the ORF has to do with the appearance of 
derivatives of conditional probabilities, which for deterministic systems would correspond to 
groups of very close trajectories. These first appear in iS*-^-*, see for instance Eqs. ( l45l) for 
T^*^^-* . The nonlinear character of the KLRF means that such terms appear in comparatively 
higher order of the perturbation theory. For example, iS^^-* contributes to T?,^^-* but not to 

(3) 

the Fisher information. Instead it first contributes to Qi23- 

We have demonstrated that the Fisher information behaves in a qualitatively different 
way depending on whether the system is perturbed from equilibrium or from an out-of- 
equihbrium steady state. For the classically stochastic systems considered here we have 

('2) 

seen that in the former case Q12 (^Ij'^s) = f{ti + 2^2)- That is, the Fisher information is a 
one dimensional function of the two time delays. This qualitative difference results from the 
self-duality of equilibrium dynamics, which is another expression of the principle of detailed 
balance. Q^^a does not show a similar reduction of dimension [see Eq. (IB6|1 ]. This property 
is special to the Fisher information. 

We have focused on the properties of the Fisher information for overdamped stochastic 
dynamics. Do other types of systems exhibit similar behavior? In App. |A] we consider 
deterministic Hamiltonian systems. In that case the unitary dynamics results in a Fisher 
information which depends only on the time ti. It will be of interest to study how (not 
overdamped) stochastic systems bridge between the overdamped and unitary limits. 

We showed that the KLRF can serve as a useful measure characterizing the system's 
dynamics. They encode information which differs from the information encoded in the 
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ORF. This is demonstrated by the abihty of the Fisher information to distinguish between 
systems perturbed out of equihbrium or out of a non-equihbrium steady state. We expect 
other useful properties of the KLRF to be revealed by further studies. 
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Appendix A: Unitary dynamics 

For completeness, in this Appendix wc discuss the application of nonlinear response 
theory to deterministic Hamiltonian systems. We use simple examples to clarify the relation 
between the response of a system and the KuUback-Leibler distance. The Hamiltonian of 
the system is assumed to be of the form H — Hq — s{t)A, where A is the perturbation. We 
start with some general comments stemming from the fact that the dynamics in phase space 
is an incompressible flow. 

1. General considerations 

The state of a classical system is described by its phase space density p. It is important 
to note that this is the full probability density, which includes information on all the degrees 
of freedom. Let us denote the propagator of the classical trajectories by h, so that x(t) = 
/It-(x(0)) where x denotes a phase space point (all coordinates and momenta). Similarly, 
the propagator for the probability distribution is denoted by W(r). Liouville theorem tells 
us that phase space volumes do not change in time. As a result, there is a simple relation 
between the phase space density at different times, 

p(x, r) = W(r)p(x, 0) = p(/i-.(x), 0) = p(x(0), 0). (Al) 

The density is simply transported with the dynamics in phase space. The density at hr^x) 
at time r is equal to the density at x at time 0. 
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This property of the phase space dynamics have an interesting consequence. Any integral 
whose integrand depends locally on p alone is time independent, since the values of the 
integrand are just transported around by the dynamics. An interesting example is the 
entropy function 

5(r) = -y"rfxp(x,r)lnp(x,r). (A2) 

Let us change the integration variable to xq = /;,_r(x), which is just the phase space point 
which would flow to x after a time r. Liouville theorem assures us that the Jacobian of the 
transformation is unity and therefore 

S{t) = - j d^oPi^o, 0) lnp(xo, 0) = So, (A3) 

where we have used Eq. flAl|) . It is clear that this entropy does not depend on time. 

Eq. (IA3I) is a result of the unitary evolution in phase space. It connects with all the 
intricate problems related to the emergence of macroscopic irreversibility out of microscopic 
reversible dynamics. This deep problem is beyond the scope of the current paper. Such 
problems are circumvented when one uses a reduced probability density, whose dynamics is 
irreversible to begin with. 

Unitary dynamics lead to an interesting result for the reverse KuUback-Leibler distance 

X'(pIIpq) = Tr pin — = Tr plnp — Tr pin po. (A4) 

Po 

We have seen that the first term is a constant of the motion. Assuming a Hamiltonian 
system intially in equilibrium po = e~^^"/Z. We have 

V{p\ Ipo) = -So + Tr pf3 {Ho - J'o) = P Tr (p - po) -Ho, (A5) 

where J-q = —ksT In Z is the free energy of the initial equilibrium state, '^'(pl |po) is therefore 
linear in p, and thus it is equivalent to a calculation of response functions! 

The general comments above raise two interesting points. First, the fact that terms such 
as Tr p In p are constant means that their expansion in powers of pulse strengths Si turns 
out to have only the constant term, all other terms in the expansion must vanish. 

The second point of interest has to do with the relation between the Kullback-Leibler 
distance and its reverse. Generally, V{p\\po) ^ T){po\\p\ We are interested in systems 
pertubed by a series of pulses and compare the initial and final distributions. In that case 
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we can use unitarity to change variables from the phase space points at the final time to the 
points at the initial time which are connected to it by the dynamics. This gives 

V{p\\po)= [ dM^,T)\nP^= [ ^;xopo(xo)ln4^ = I?(po||p). (A6) 
J Po(x) J p(xo,0) 

Here Po(x(t)) = p(xo,0), that is, p is the density that would evolve to the equihbrium 

density under the influence of the pulses. It is not equal to p(x, r), which tells us that the 

distance and its reverse are not equal. 

While in general the distance and its reverse are not equal, when the distance between 

the distributions is small, its easy to show that their leading order expansion in 5p = p — po 

is the same. Loosely speaking 

I?(po||p) = IJ d^(^^ + 0{6p')^ (A7) 

^^(pIIpo) = IJ rfx + 0(5p3)) . (A8) 

(A9) 

We can deduce that the Fisher information, which is determined by this leading order, could 
be obtained from a calculation of response for systems with unitary evolution. (Since it could 
also be obtained from the reverse distance.) 

Our last general point is also related to unitarity, but has to do with the Fisher infor- 
mation. According to Eq. f l22|) the Fisher information is built from two partial densities. 
(More accurately, density differences.) These evolve with respect to the same Hamiltonian 
"Ho between interactions with the pulses. As a result, the integral in Eq. fl2^ do not depend 
on t2, for the same reason that caused the entropy S to be time independent. We flnd that 

QS(t2,ti) = Ql?(ti), (Aio) 

as a result of the unitarity of the phase space dynamics. Similar independence on the flnal 
time interval would also appear for higher order terms in the expansion of the Kullback- 
Leibler distance. 



2. A single Harmonic oscillator 

The simplest system that can serve as an example is a single Harmonic oscillator 

^0 = ^ + ImQ'Q'. (All) 



21 



We take the perturbation to be 

A{Q) = ae-^l^\ (A12) 
For this system it is trivial to solve for the free evolution 

g(r) = g(0)cos^]r + :^sinfir 

P(r) = P(0)cos^]r-^]Mg(0)sin^]r. (A13) 

This relation can be inverted, expressing P(0),g(0) in terms of P{t),Q{t) 

P(t) 

Q(0) = Q(r) cos ^]r ^sin^]r 

Mil 

P(0) = fiMg(r) sin ^]r + P(r) cos fir. (A14) 

Eq. (IA14P will be useful when one propagates probability distributions in time. 
The equilibrium distribution of the Harmonic oscillator is a Gaussian 



Po = — ^ — e 2M 2P^>"' ^ . (A15) 



We also note that 



Capo = {A, po} = T7^^e-«/«Vo, (A16) 

is the linear order correction for the density just after interaction with one pulse. 

To calculate S^^^ we need to propagate Eq. ( lAlGll . With the help of Eqs. (lAll) and ( ]A14p 
we find 

S^^\Q,P;ti) = -^(fiMgsinfiti + Pcosfiti)e~^['^™'''*'~^""™']po(Q,i'). (A17) 
MQq 

In the derivation we have used the fact that po is invariant under the evolution with respect 
to T-Lq. 

To calculate 5*^^^ we operate with Ca on S^^\ and then propagate the resulting correction 
for the density for a time interval t2. A straightforward calculation leads to 



5(2)(g,P;t2,ti) = Uih)CAS^'^ 



MQl 



cos Vtti + {VtMQ sin VL{ti +t2) + P cos fi(ti + ta)) 



-— [nMQ sin nt2 + P cos nt2] + — ^—sinfiti ) I (A18) 

g p 

X exp -— (cos Qt2 + cos Q{ti + ta)) + „ „ (sin Qt2 + sin r2(ti + ta)) Po 

L vo ikzizgo J 
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We now turn to calculate the first two ORF, using Eqs. (12T!) and (123|) . The calculation 
of TZ^^^ is cumbersome but straightforward. 



dQdPA{Q)S^^\Q,P;t 



a 



QMQi 



sin Qti exp 



^1 + cosilti' 



(A19) 



The calculation of 7^'-^^ is more involved, and we only include the final result 



a 



X exp 



sin ^2^2 [sin Qti + sin Q{ti + t2)] 



1 



h cos fiti + COS VLt2 + COS VL{ti + t2 



. (A20) 



We would like to compare these response functions to the KLRF, and in particular to 
the Fisher information, which can be calculated using Eq. f l22|) . The calculation can be 
simplified by using the classical coordinates right after the second interaction with the pulse 
as integration variables. Due to unitarity, the Fisher information only depends on the time 
difference ti, see also the discussion in the previous subsection. 

We get 



QS?(t2,ti) 



cosr2ti — 



sin^ Vtti 



exp 



(A21) 



This expression seems similar to the first order response function. One can indeed show that 
they are related by 



d 



(A22) 



It will be interesting to check whether this expression could be generalized to any Hamilto- 
nian system (with unitary dynamics). 

Appendix B: Third order ORF and KLRF 

In the main text we have calculated the ORF and KLRF to second order. In the appendix 



we present the third order quantities, 7l^^\ and Q 



(3) 
123- 



1. Perturbation theory 

Following the calculations performed in Sec. [Tll the third order response functions are 
j;iven by 

7^^'^(^3,^2,^l) = I dxA{x)S^^\x-hMM), (Bl) 
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and 

+5(1) (x; h + h)S^^^ {x- h, h + h) + 5(1) (x; h)S^''^ (x; ^3 + ^2, h)] 

-2pf{x)S^^\x- ti+h + h)S^^\x- h + h)S^^\x- h)] , (B2) 

where 

5(^)(x;t3,t2,^l) = y rfXi ■ ■ ■ dXQUo{x,X^;t'i)CA{Xfi,X^)UQ{x^,Xi,t2)CA{Xi,X'i) 

X Wo(a;3, a;2; ti)£A(a;2, xi)pq{xi). (B3) 

2. Eigenvalue expansions 

Here we present expressions for the third order ORF and KLRF in term of the eigenvalues 
and eigenfunctions of the stochastic dynamics. At high orders the nonhnear character of the 
contributions to the Kullback-Leibler distance result in integrals with products of several 
eigenfunctions. Here we will only need the one with three eigenfunctions 

Jnmi = j dxqn{x)qm{x)qi{x)po{x). (B4) 

The third order ORF is given by 

7^('^^3,t2,tl) = 5^C^^i3,„S„.„i3„oe-^(")*^-^('")*^-^(')*^ (B5) 

Imn 

This ORF should be compared to the third order KLRF, which is calculated using Eq. flB2p . 
We find 

QfkhMM) = Y^^^^^^uBno {e-^(")*^e-'^(-)[*^+*^+2*3] 

mn 

_^g-A(n)[ti+t2]g-A(m)[t2+2t3] _^ g-A(n)ti g-A(m)[t2+2t3] | 

- 2 5^i3„oS^oi3,oJnwe-^(")[*^+*^+*^le-^('")[*^+*^]e-^(')*^ (B6) 

nml 

The expression for presented in Eq. (IB6p . has three terms in which the orthogonality 

condition (!5T|) has been used, pointing to a reduction of dimension in the time dependence 
of this specific term. However, the time combinations in these terms are all different. In 
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addition, the fourth term in Eq. f lB6|) clearly depends on all its time variables. We conclude 
that in contract to the Fisher information, the higher order KLRF depend on all their time 
variables. The reduction of dimension is therefore specific for the Fisher information. It 
results from the fact that it is built out of a single product of two density corrections. 



Appendix C: The overdamped harmonic oscillator 



In this Appendix we consider a simple example of an overdamped harmonic oscillator 
with a potential l/2Mf2(5^, with a perturbing potential A{Q) = OLe~^l'^° . For this system it 
is possible to write explicit expressions for the eigenvalues and eigenfunctions of the Fokker- 
Planck operator, as well as to perform several of the integral, defined in Sec. Illli In this 

case 



. 1 d 
^oP = 7^ 



8 ' 

MQ'^Qp + keT—p 



dQ' 



Cap 



a d 
Qol dQ 



(CI) 
(C2) 



This model has been studied extensively 22|]. The equilibrium density is given by 



Po(Q) 
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The left eigenfunctions, and the eigenvalues, of this model are 



(C3) 



qn{Q) 



2"n! 



A(?7,) = n. 

7 

Here, Hn are the Hermit polynomials. 

The following properties of the Hermit polynomials 



(C4) 
(C5) 



dx 



2nHn-i{x) 



HJx) = (-W 



dx^ 



(C6) 
(C7) 



together with the exponential form of A{Q), allow to compute some of the integrals defined 
in Eqs. ([52]) - (jBi]) explicitly. 
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For instance, 



CAn = aj rfQe-«/««g„(Q)po(Q) = «Y^/ dye-^^f^ HMe-'" ■ (C8) 

One can substitute Eq. flC7p for the Hermit polynomial, and use integration by parts. This 
leads to 



n/2 

Cau = a{-lT\l--Q-,- ( ^ ) / dye-y^-^^^"^. (C9) 



This is a Gaussian integral which is easily evaluated to give 

n/2 fe 



One can also easily calculate the integrals 



S„„ . _ /dQ|^|ip„W) ^ ^ [,yH..M^"%-*'/^, (Cll) 



by the same technique. One finds 



Bno = (-1)"-^^-^ ( ^ V' (C12) 



Comparing Eqs. (IClOP and (]C12p . we see that for this model 



Bno = -n——^CAn. (C13) 



As a result, there is a simple relation between the Fisher information and the first order 
response function, 

(C14) 



This relation is a special result for this model, and is not expected to hold for other systems. 

One can also obtain explicit results for Bnm with non vanishing indices. However, the 
calculation and the result are quite combersome, and are omitted. We were not able to 
calculate Jnmx explicitly. 
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